Damped simple pendulum

In this notebook we will solve the equation of motion for a simple pendulum with damping.

Libraries

We will use the following libraries:

  • numpy for numerical calculations
  • matplotlib for plotting
  • scipy.integrate for numerical integration
  • ipywidgets for interactive plots
  • IPython.display for displaying the plots
  • sympy for symbolic calculations
In [1]:

Mathematical model

For this system, we will consider a simple pendulum with damping. So we have the following parameters:

  • m is the mass of the pendulum bob
  • l is the length of the pendulum
  • g is the acceleration of gravity
  • b is the damping coefficient
  • θ is the angle of the pendulum
  • ω is the angular velocity of the pendulum
  • t is the time
In [2]:

Now we need to define the position x and y of the pendulum bob. We can do this by using the following equations:

x=lsin(θ)y=lcos(θ)

after that, we can define the kinetic energy T and the potential energy V of the system:

T=12m(x˙2+y˙2)V=mgy

In [4]:
Out[4]:
(0.5m(l2sin2(θ(t))(ddtθ(t))2+l2cos2(θ(t))(ddtθ(t))2), glmcos(θ(t)))

Now we can define the Lagrangian L of the system:

L=TV

In [5]:
Out[5]:
lm(gcos(θ(t))+0.5l(ddtθ(t))2)

As we get the Lagrangian, we can use the Euler-Lagrange equation to get the equation of motion, but we use them for a system with non-conservative forces, so we need to use the Euler-Lagrange equation for non-conservative forces:

ddt(Lq˙j)LqjQj=0

Qj=i=1NFiriqj

where qj is the j-th generalized coordinate, q˙j is the time derivative of the j-th generalized coordinate, ri is the position vector of the i-th particle, Fi is the force vector acting on the i-th particle, and N is the number of particles.

In [6]:
Out[6]:
bl2ddtθ(t)+glmsin(θ(t))+l2md2dt2θ(t)

For simplicty we can divide the equation by l2m

In [7]:
Out[7]:
bddtθ(t)m+gsin(θ(t))l+d2dt2θ(t)

Now that we have the equation of motion, we can solve it numerically.

In [8]:

Let's define a vector S that returns dSdt.

In [11]:

Now we can solve the equation of motion numerically.

In [12]:
In [13]:
0510152025303540-2.5-2.0-1.5-1.0-0.50.00.51.01.52.02.53.0Time (s)Angle (rad)Angle vs Time0510152025303540-6-5-4-3-2-10123456Time (s)Angular velocity (rad/s)Angular velocity vs Time

We can plot the phase space of the system.

In [14]:
-2.5-2.0-1.5-1.0-0.50.00.51.01.52.02.53.0-6-5-4-3-2-10123456Angle (rad)Angular velocity (rad/s)Phase space

If we want to see a plot of the phase space for diferent initial conditions, we can use streamplot.

In [15]:
-20-15-10-505101520-20-15-10-505101520Angle (rad)Angular velocity (rad/s)Phase space

Finally, we can plot the position of the pendulum bob.

In [16]:
Out[16]:
-1.4-1.2-1.0-0.8-0.6-0.4-0.20.00.20.40.60.81.01.21.4-1.4-1.2-1.0-0.8-0.6-0.4-0.20.00.20.40.60.81.01.21.4time = 40.0sx (m)y (m)Pendulum
In [ ]: